Vielleicht ist es möglich Objekte aus dem IR-10.8-µm-Kanal über die Minima und davon ausgehenden Grenzen zu definieren.
Um das zu untersuchen, brauchen wir zunächst ein paar Pakete.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage as ndi
from skimage.feature import peak_local_max
import cv2
import datetime as dt
import l15_msevi.msevi as msv
import sys
sys.path.append("/vols/talos/home/stephan/utils/tracking")
import tracking_common as tc
Wir laden für den 28.07.2012, 12:00 UTC ein Beispiel.
t = dt.datetime(2012,7,28,12,0)
s = msv.MSevi(time=t,chan_list=["IR_108"])
s.lonlat()
s.rad2bt()
ir108 = s.bt['IR_108']
Nun suchen wir nach den lokalen Minima im Bild.
minima = peak_local_max(-ir108, min_distance=2,threshold_abs=-273)
print minima
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ax.imshow(ir108,vmin=210,vmax=300,cmap='gray_r')
x = [m[1] for m in minima]
y = [m[0] for m in minima]
ax.scatter(x,y,marker='+',color='blue')
Das sind ziemlich viele Minima. Es gibt mehrere Möglichkeiten deren Anzahl zu reduzieren:
Das probieren wir im Folgenden mal aus.
Für die Glättung gibt es mehrere Möglichkeiten. Hier probieren wir Folgende aus:
All diese Verfahren sind abhängig von der Filterbreite. Was ist ein guter Wert dafür?
ir108_transformed = (tc.scale_array_min_max(ir108)*255).astype("uint8")
mean_blur = cv2.blur(ir108_transformed,(5,5))
gauss_blur = cv2.GaussianBlur(ir108_transformed,(5,5),0)
median_blur = cv2.medianBlur(ir108_transformed,5)
bilateral_blur = cv2.bilateralFilter(ir108_transformed,15,15,15)
plt.imshow(ir108_transformed)
plt.colorbar()
fig,ax = plt.subplots(2,2,figsize=(14,14))
ax[0,0].imshow(mean_blur,vmin=0,vmax=255,cmap='gray_r')
ax[0,0].set_title(u"Mittelwertglättung")
ax[0,1].imshow(gauss_blur,vmin=0,vmax=255,cmap='gray_r')
ax[0,1].set_title(u"Gaußglättung")
ax[1,0].imshow(median_blur,vmin=0,vmax=255,cmap='gray_r')
ax[1,0].set_title(u"Medianglättung")
ax[1,1].imshow(bilateral_blur,vmin=0,vmax=255,cmap='gray_r')
ax[1,1].set_title(u"bilateraler Filter")
Die Glättung mit dem bilateralen Filter sieht subjektiv am besten aus. Das Feld ist relativ stark geglättet, aber die Kanten sind noch gut erhalten.
thresh = (273 - np.min(ir108)) / (np.max(ir108) - np.min(ir108))*255
thresh
minima = peak_local_max(-bilateral_blur, min_distance=2,threshold_abs=thresh)
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ax.imshow(ir108,vmin=210,vmax=300,cmap='gray_r')
x = [m[1] for m in minima]
y = [m[0] for m in minima]
ax.scatter(x,y,marker='+',color='blue')
Das sieht schon besser aus. Nur die Häufungen von Minima sind seltsam.
minima = peak_local_max(-ir108, min_distance=5,threshold_abs=-273)
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ax.imshow(ir108,vmin=210,vmax=300,cmap='gray_r')
x = [m[1] for m in minima]
y = [m[0] for m in minima]
ax.scatter(x,y,marker='+',color='blue')
Das sieht besser aus.
minima = peak_local_max(-ir108, min_distance=2,threshold_abs=-250)
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ax.imshow(ir108,vmin=210,vmax=300,cmap='gray_r')
x = [m[1] for m in minima]
y = [m[0] for m in minima]
ax.scatter(x,y,marker='+',color='blue')
Das sind immernoch viele Minima, aber sinnvoller verteilt.
th = 250
thresh = (th - np.min(ir108)) / (np.max(ir108) - np.min(ir108))*255
test = 255 - bilateral_blur
thr = 255 - thresh
minima = peak_local_max(test, min_distance=15,threshold_abs=thr)
#from skimage.morphology import extrema
#minima = extrema.h_minima(bilateral_blur,30)
minima
test = 255 - bilateral_blur
from mpl_toolkits.axes_grid1 import make_axes_locatable
def colourbar(mappable):
ax = mappable.axes
fig = ax.figure
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
return fig.colorbar(mappable, cax=cax)
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
#ax.imshow(np.ma.masked_less(minima,1))
test_plot = ax.imshow(test,vmin=0,vmax=255,cmap='gray')
colourbar(test_plot)
print(np.max(ir108) - np.min(ir108))
print(th - np.min(ir108))
print thresh
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ax.imshow(ir108,vmin=210,vmax=300,cmap='gray_r')
x = [m[1] for m in minima]
y = [m[0] for m in minima]
ax.scatter(x,y,marker='+',color='blue')
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
test_plot = ax.imshow(np.ma.masked_less(test,thr))
colourbar(test_plot)
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
test_plot = ax.imshow(np.ma.masked_greater(ir108,250))
colourbar(test_plot)
Das sieht schon recht gut aus. Bloß dies Häufung von Minima ist seltsam. Das kann man wahrschenlich folgendermaßen beseitigen:
from scipy import ndimage as ndi
new_minima = np.zeros_like(ir108)
for m in minima:
new_minima[m[0],m[1]] = 1
lab,nlab = ndi.measurements.label(new_minima)
new_mins = ndi.center_of_mass(new_minima,lab,np.arange(1,nlab))
new_mins
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ir_plot = ax.imshow(ir108,vmin=210,vmax=300,cmap='gray_r')
colourbar(ir_plot)
x = [m[1] for m in new_mins]
y = [m[0] for m in new_mins]
ax.scatter(x,y,marker='+',color='blue')
Das sieht als repräsentativer für die Objekte aund Startpunkte für die Objektdefinition ziemlich gut aus. Es sind Minima von nichtkonvektiven Wolken dabei, aber um diese zu entfernen, ist wahrscheinlich eine Wolkentypisierung nötig.
ct = msv.MSGtools.get_nwcsaf_prod('CT',t,calibrate=True)
ct = np.ma.masked_where(np.logical_or(ct<5, ct>10),ct)
h = 14
fig,ax = plt.subplots(1,1,figsize=(h*1.61803,h))
ct_plot = ax.imshow(ct)
colourbar(ct_plot)